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Abstract. The interaction force between likely charged particles/surfaces is usually 
repulsive due to the Coulomb interaction. However, the coimterintuitive like-charge 
attraction in electrolytes has been frequently observed in experiments, which has been 
theoretically debated for a long time. It is widely known that the mean field Poisson- 
Boltzmarm theory cannot explain or predict this anomalous feature since it ignores 
many-body properties. In this paper, we develop efficient algorithm and perform the 
force calculation between two interfaces using a set of self-consistent equations which 
properly takes into accoimt the electrostatic correlation and the dielectric-boundary ef¬ 
fects. By solving the equations and calculating the pressure with the Debye-charging 
process, we show that the self-consistent equations could be used to study the attrac¬ 
tion between like-charge surfaces from weak-coupling to mediate-coupling regime, 
and that the attraction is due to the electrostatics-driven entropic force which is signif¬ 
icantly enhanced by the dielectric depletion of mobile ions. A systematic investigation 
shows that the interaction forces can be tuned by material permittivity, ionic size and 
valence, and salt concentration, and that the like-charge attraction exists only for spe¬ 
cific regime of these parameters. 

PACS: 82.45.Un, 64.70.pv, 82.60.Lf 
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1 Introduction 

The anomalous attraction between likely charged particles has been widely discussed 
during the past decades since the experimental observation in colloidal BT^ and 
biological systems |1U|23. The like-charge attraction (LCA) at strong-coupling regime 
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(high surface charge, multivalent counterions or low temperature) has been reported by 
particle-based simulations and theoretical studies such as density functional theory and 
strong-coupling theory Ill8ll32ll35ll39ll43tl4^ . Colloidal suspensions in monovalent elec¬ 
trolytes are usually considered as systems in weak-coupling regime at which the mean- 
field Poisson-Boltzmann (PB) theory is generally useful. The forces between charged 
surfaces are usually modelled through the Derjaguin-Landau-Verwey-Overbeek (DLVO) 
theory, which combines the linearized PB equation with the Van der Waals interaction. 
However, the PB theory always predicts repulsive interaction between like charges Il48l 
l56ll . and the DLVO theory overestimates the repulsion or underestimates the attraction 
Il28ll . in contrast to the results of LCA experiments for colloids IIMIISSH . which demon¬ 
strates many-body effects are still significantly important at weak-coupling conditions 
and can not be ignored in order to understand the LCA phenomenon systematically. 

For two colloidal particles with a small separation, the depletion force induced by the 
entropic repulsion of smaller mobile ions in the solvent could result in an effective inter¬ 
colloid attraction. The excluded volumes of small particles lead to a depletion layer. 
When the separation between surfaces is at the range of the particle diameter, they are 
attractive due to this depletion interaction ||2[|55l. In electrolytes, this attractive distance 
can be much larger, comparable with the Debye screening length due to the long-range 
nature of electrostatic interaction IITSlI . In this sense, the attractive Van der Waals and 
the Casimir force HTTHsTIISTII can be ignored in the LCA analysis. In the presence of in¬ 
homogeneous dielectric permittivity (whose properties have attracted wide recent inter¬ 
est II51I16II17II20II30II59II63II L ions are repelled from the low dielectric surfaces by their image 
charges, which has been shown in particle-based simulations II111I24II . When the separa¬ 
tion of surfaces becomes narrower, the repulsion becomes stronger, resulting in lower 
ionic concentration between the surfaces, and thus leading to an entropic driven attrac¬ 
tion. The electrical field can also influence the solvent alignment and dielectric property, 
which is considered to be important for high ionic concentration or at the strong coupling 
region. This solvent polarization effect could be modelled through the Langevin theory 
to yield a field-dependent dielectric coefficient BllTlISlI . Recently, a continuum theory with 
a set of self-consistent equations including fluctuation effects of ions is introduced as a 
variational approach with general Gaussian ansatz II47[I58II . which have been generalized 
to take into account the image charge effects on attractive forces between neutral plates 
and in weak coupling region II57II . In the theory, the self energy of mobile ions is used as a 
correction to the mean potential in the PB theory, leading to a more accurate approxima¬ 
tion of the potential of mean force (PMF) in the Boltzmann distribution, which consists 
of contributions from both the local ionic correlation and the ion-interface interaction en¬ 
ergy. Although the LCA phenomenon has been predicted, it is less understood what is 
the effects of salt property as well as the charged surface. 

In this paper, we study the interaction between plates in electrolytes by the self- 
consistent field model developed in IlHlIfflU which includes a modified treatment of the 
self energy. Through the modification, we include the ionic size effect to avoid the ion 
collapse in the original theory and are able to study systems from weak-coupling to 
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mediate-coupling regimes and with a inhomogeneous dielectric coefficient beyond the 
point-charge models. We present a systematic analysis for effects of the surface charge 
density, the bulk salt density, and salt species on the depletion-induced LCA. It is shown 
that for symmetric electrolytes the interaction force between plates is repulsive when the 
surface charge density is high or when the salt concentration is low. The force becomes 
attractive at the opposite regimes, for which the strength will be stronger when the charge 
density is decreased or the salt density is increased. 


2 Self-consistent field model 


Consider an equilibrium charged system with N species of mobile ions in solution. Given 
the PMF of the zth species, lUi (r), the ion distribution can be determined by the Boltzmann 
distribution, 

c;(r) = cfe-^“'-«, (2.1) 

where is the bulk concentration, and j6 = 1 / kg T is the inverse thermal energy with 
ks the Boltzmann constant and T the absolute temperature. The PMF describes the free 
energy change by moving the ion from the bulk region into the current position, which is 
often given by particle simulations or advanced theory in statistical physics. If the PMF 
is approximated by the mean potential energy W; = z,e^, with ^ being determined by the 
Poisson equation, 

N 

-eoV-eVO=^z,-ec,, (2.2) 

! = 1 


where Cq is the vacuum dielectric constant, e is the relative permittivity of the media, and 
Zi is the valence, we obtain the PB equation since C; is now an explicit function of O. 

To properly take into account the dielectric polarization and the correlation effects 
of mobile ions, one should study the self energy or the intrinsic chemical potential of 
a test ion, which is considered as a correction term to the mean potential to better ap¬ 
proximate the PMF |l3l[T3[T3H29lllQlll3llZll5Ql[57K^. In our model, we treat an ion as an 
ion-inaccessible sphere of radius fl/ with the point charge z, at the center. This allows us to 
deal with variable dielectric media where the ionic Born energy, which strongly depends 
on the ion radius, is not constant and thus cannot be discarded. For a homogeneous sys¬ 
tem, the self energy of an ion can be described by the exact solution of the Debye-Hiickel 
equation IIMlI . 


u,=^ui= y \ 

2 Sneoel+KUi 


(2.3) 


where k is the inverse Debye length of the electrolyte. In the vicinity of a surface, the 
ionic distribution and the dielectric permittivity are space-dependent, and a direct use of 
Eq. (|2.3|) by replacing k by the space-dependent quantity cannot account for the nonlocal 
property of the electrostatic correlation. We assume that the ion size is small and that the 
presence of a test ion does not influence the ionic distribution in the electrolyte except for 
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the region occupied by the ion due to the finite excluded volume. The self energy of the 
ith ion species is then formulated through the generalized Debye-Hiickel equation 1^ . 


-£oV-£,VG/(r,r')+2I,G,(r,r')=^(r,r'), 
u,(r) = lim[G,(r,r') -Go(r,r')], 

r'—>r 


(2.4) 


where the dielectric permittivity and the ionic strength locally depend on the site of the 
test ion. 


with 


l£(r), otherwise. 

(2.5) 

, fo, if |r-r'|<fl„ 

T(r,r ) = < 

1 f(r), otherwise. 

(2.6) 


(2.7) 


being the mean ionic strength. The function Go in Eq. (12.41) is the Green's function in free 
space, satisfying 

-£oV2Go = (i(r,r'), (2.8) 


which is used to eliminate the singularity in G/. 

When £ and I are constant, the solution iq reduces to the Debye-Hiickel theory (12.31) . 
The equation for G;, Eq. (I2.4I) . being considered as the generalization of the Debye- 
Hiickel equation is in the sense that the inverse Debye length depends on the space, 
K:(r) = y^f(r)/27r£o£w/ where £w is relative dielectric constant of water. With the expres¬ 
sion for the self energy, we then obtain a modified PB equation. 


N 

— £oV •£VO= Y^ZiCc’-exp 







(2.9) 


where exp [jiz^e^u^ /2) and w■’ is the bulk value of m,. 

So far, we have described a continuum model for inhomogeneous electrolytes through 
a set of self-consistent field (SCE) Eq. (I2.4I) - (I2.9D . This set of equations include the long- 
range correlation by the solution of the generalized Debye-Hiickel equation. However, 
the model may make imprecise prediction for dense electrolytes or strongly-correlated 
systems, because the hard-core correlation is neglected. In spite of this, it was reported 
that the SCE model can capture correlation phenomena well for electrolyte systems up to 
the coupling parameter H~50 ll4T]| . which is certainly in the regime of strong coupling. 
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Figure 1: A schematic description of electrolytes between two likely charged planes of separation D. Ions are 
of finite sizes and are inaccessible to a zone of thickness from the surfaces. If the dielectric constant outside 
the electrolyte is much smaller than that of the solvent, a repulsive self energy will be acted on each mobile ion. 


3 Method 

3.1 Physical setup 


We compute the interaction forces between two charged surfaces in a system schemati¬ 
cally shown in Fig. [ij The two parallel planes are with surface charge densities and 
(72, respectively. The planes are orthogonal to the x axis, locating at x = ztDfl. An elec¬ 
trolyte is in between the two planes, and its relative dielectric permittivity is ew The 
outside space is filled with dielectric media characterized by the relative dielectric con¬ 
stant Sb- We assume negative surface charges on two planes, and study electrolytes with 
two species of ions (counterions are positive ions). Since mobile ions have finite sizes, 
there is a zone, i.e. the Stern layer, of the ion-radius thickness inaccessible by mobile ions 
near each charged plane. 

Due to the planar geometries, the ionic density is invariant along the y — z plane par¬ 
allel to the surfaces. The modified PB equation (I2.9D is one-dimensional, and can be it¬ 
eratively solved by a finite-difference method, giving the self energy. In the calculation, 
Neumann boundary conditions ^ = 0 are applied on two ends of interval [—L,L] with 
L > D / 2 big enough to remove artifacts from the boundary. The surface charge placed at 
±D/2 is approximated by the Kronecker delta with a coefficient inversely proportional 
to the mesh size. 











6 


3.2 Numerical method 

The generalized Debye-Hiickel equation (12.41) is high dimensional and hard to solve. In 
order to develop an efficient numerical method, we first employ an asymptotic treatment 
EH, which approximately separates the self energy as the sum of two contributions, 

= + (3.1) 


under the assumptions that the ionic radius fl; can be considered as a small parameter and 
that the ionic strength is slowly varying. The first term in Eq. (|3.1|) is the Born solvation 
energy Il9ll58ll, 


,Born 


47reo«! 


1-1 
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(3.2) 


The second term containing the contribution from the long-range Coulomb correlation is 
expressed as. 


.Cold , 


l + fl;Kf,' 


(3.3) 


where K], is the inverse Debye length in the bulk. The value of u is determined by the 
following equations. 


I -eoV-eVG(r,r')+21G(r,r') = ^(r,r'), 
lM(r) = lim[G(r,r')-Go(r,r')]. 

r'^-r 

This formulation (13.31) is similar to the description of the intrinsic chemical potential in 
Ref. II53II . The form of —K},/{1+Ki,a) as the correlation energy was also established by 
the means of integral equation theory EH- For homogeneous system 1 = Kj£oew/2 and 
e(r) =ew/ it is easy to see u = —K;,/47r£oew- However, for systems with strong inhomo¬ 
geneity, such as those in our results, u can be largely different from —Kjj/Aneoe. For such 
systems, we find our asymptotic expression (13.31) works well. In the previous work BTH . 
another expression ^u/{l — 4iTzeoeaju) was used, which also works if u does not 
become positive. In this study, u could become a large positive number near charged sur¬ 
faces, thus this is essential to use expression (13.3b in order to avoid the zero-denominator 
instability. 

Eq. (13.4b can be derived through field theoretical approach and has been studied as 
the self energy of a point charge. The advantage of Eq. (I3.4|) compared to Eq. (I2.4|) is 
that it removes the multiple length scales in the model and can be numerically computed 
efficiently MMlI . Then the equation can be better understood if we write it into the form: 
AG=5 where A = —eoV-eV-|-2f is an operator which can be viewed as a matrix and 5 can 
be viewed as the identity matrix. Since we are only interested in the diagonal entries of G, 
the selected inversion algorithm II37II could be used and it is very efficient. Considering 
the geometric symmetry, we employ the Fourier transform in the y—z coordinates to the 
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Green's function equation. 


3 3 2 r^Tl \ 

-eoT-e—+eo£CL>+2I{x) 
dx dx 


G{u};x,x') = —5{x-x'). 


(3.5) 


In domain x G [—L,L], we could discrete the equation by finite-difference method with 
periodic boundary conditions, which is in the form of: A„G = where A„ is a nxn 
symmetric sparse matrix and I„ is the identity matrix. Then the diagonal part of G can 
be solved by the selected inversion method, whose computational cost is at the order of 
Cholesky factorization. The divergent part of the self-Green's function is eliminated by 
numerically solving Go{co}x,x) by the same scheme, then G —Go is convergent with the 
refinement of the mesh size. The self energy is then determined by the inverse Fourier 
transform ISTH . 

We should also notice that, as Go ~ 1/ |r—r'| decays very slowly, we can hardly guar¬ 
antee the precision of Go without choosing relatively large computational domain. To 
avoid this inconvenience, instead of solving —eoV^Go = ^i(r,r'), we solve — V^Go+k^Go = 
<5(r,r')/£o where A:>0 is a constant and use the fact that limj/^r(Go — Go) = —k/Anso- 

Now we have described the method for solving the self energy of ions with given 
ionic distributions. With the numerically computed self energy, Eq. (12.9b for the electric 
potential is discretized through the finite difference method of second order accuracy. 
Since the ionic distributions are nonlinear functions of the electric potential, we employ 
the quasi-Newton method for the iteration in solving the nonlinear PDFs. Then the up¬ 
dated ionic concentrations are given by Boltzmann's distribution and the whole system 
is solved through the fixed-point iteration until the final self energy and electrical po¬ 
tential are self-consistent to the required precision (error less than 10^^). After solving 
the set of self-consistent equations, we obtain the equilibrium electrical potential, ionic 
concentrations as well as the self energy. 


3.3 Free energy and pressure 


In order to measure the interaction force between the two planes, we first calculate the 
grand potential per unit area as the function of the separation. 


f(D) = 


rD/2 

-D/2 


525%i£+WEc,(log|-l)+/n 


dx, 


where /a is the contribution of fluctuation, given by the Debye charging process 


ffi = 'Ezfe^Ci / A 


Ua 


l+ciiXKi, AneoUi \£ 


■-1 


dA, 


(3.6) 


(3.7) 


and Ma is determined by, 

r -£oV-eVGA(r,r')+2A2lGA(r,r')=d(r,r'), 

1 WA(r)^lim[GA(r,r')-Go(r,r')]. 
t r'H-r 


(3.8) 
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The fist two terms in Eq. (I3.6|) represent the electrostatic and ideal-gas entropic contri¬ 
butions on the free energy. The third term is the correlation free energy. It should be 
emphasized that the entropy of mobile ions should have additional contribution due to 
the finite size, e.g., by the modified fundamental measure theory IISTIIS^I^ . For elec¬ 
trolytes without dense ionic particles, this contribution has minor effect. 

The integral in Eq. (13.6b is calculated through the trapezoidal rule, which has a second 
order of accuracy. The main difficulty of evaluating the free energy and pressure is to 
compute the charging process (13.7b . When A is small, due to the weak screening, the 
precision of Ga requires a large computational domain. Although no explicit solution is 
valid, we could not use the same technique as for Gq. It is noticed that the leading term of 
Ga — Go is proportional to A when A is small. So, we can choose a relatively small number 
DC, inversely proportional to the computational domain L, which could guarantee that the 
accuracy of G^. If A < a, we use the linear interpolation of Go and G^ instead of solving 
(13.8b . Thus, the integration in Eq. (13.7b can be divided into two parts: (1) for small A, 
the integration from 0 to a is integrated analytically ; (2) the integration from a to 1 is 
computed with Gauss quadrature. The error from both parts could be well controlled. 

With the grand potential, the osmotic pressure per unit area is determined by 

P = -aDF-Poo, (3.9) 

where Poo = — lim 9 dP is the bulk osmotic pressure. 

D —>-00 

4 Results and discussion 

In this section, we present numerical solutions of the self-consistent equations for the 
modified PB theory and show different aspects of the LCA phenomenon. Without spe¬ 
cial statement, we take relative dielectric constants eg = 2.5 and ew = 80. The differential 
equations are calculated in a very fine mesh (mesh size Ax=0.01nm) in order to accurately 
approximate the solutions. All systems are at room temperature, i.e., £b =0.714nm. Only 
electrolytes with two ion species are studied, we will use ± for subscript i to distinguish 
cations and anions. 

Before proceeding the calculation for the LCA, we demonstrate that the SCF model 
numerically satisfies the contact value theorem 1271 and thus the model is reliable from 
this point. The theorem for monovalent electrolytes is an exact relation between the total 
contact ionic concentration, the system pressure P and the surface charge density a, 

/(^surf) — PT , (4-1) 

where c, (xsurf) is the ionic density of the ith species on the surface at Xgurf- the classical 
PB theory satisfies this relation approximately by assuming the pressure is the ideal-gas 
pressure in the bulk. Recent study has studied the contact value theorem for the SCF 
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Figure 2: Numerically calculated contact values, (a) The profiles as function of pressure for given surface charge 
densities 0 '= [—0.01, —0.1]e/nm^. In the calculation, D increases from 0.6nm to lAnm such that P and CV are 
varied, (b) The profile as function of for given pressure. In the calculation, the surface separation is fixed to 
be 3nm so that the system pressure can be viewed as a constant equal to the electrolyte bulk pressure. 


for the point-charge case II22II . Our model is more complicated because we cannot find 
explicit expressions for both the pressure and the contact concentration. We use O.IM 
1:1 electrolytes with a± = 0.2nm and the dielectric permittivity is 80 in whole space. We 
calculate the normalized contact values numerically, CV — I^;C;(xsurf)/c,;,, as functions of 
P for given cr and for given P, respectively. The two panels in Fig. |2]show the linear 
relations for both cases, which demonstrate that the contact value theorem is fulfilled. 
Compared to the classical PB theory, the SCF model presents more accurate results since 
the approximation for the pressure with the ideal-gas pressure is not necessary. 

4.1 Symmetric electrolytes 

The pressure between surfaces depends on the difference in ionic densities between the 
bulk electrolyte and the mediated electrolyte. This is certainly relevant to the surface 
charge density a and bulk ionic number density. For symmetric electrolytes, the bulk 
ionic densities for both species are equal, c\=ct=C},. Fligher a increases the total ions 
in the mediated electrolyte, as a result the attraction between surfaces is reduced. Fig. |3] 
(a) and (b) show pressure curves as function of the effective separation for systems of 1:1 
electrolytes with uniform ionic sizes, a± = 0.2nm, where the abscissa represents the effec¬ 
tive distance between two interfaces, Dgff = D — — The surface charge density takes 
—0.01 and —0.015e/nm^, corresponding to the coupling parameter 2 = 0.032 and 0.048, 
respectively, in a weak-coupling regime. For colloid-colloid interaction, one important 

characteristic length is the Debye length defined by (eoCw^BT)- For 

a O.IM electrolyte, it is Ed = 0.96nm. From Panel (a), the pressure remains positive for 
separation much smaller or much larger than the Debye length, and shows an attrac- 
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lal {e/nm^} ^ {nm} 


Figure 3: Results for 1:1 electrolytes with uniform ionic sizes of radius a± = 0.2nm. (a) Pressure of different 
surface charge densities as function of separation distance in an electrolyte of bulk concentration c^=0.1M; 
(b)Pressure of planes with the same charge density ctj = <72 =—O.Ole/wm^ in electrolytes of different bulk 
concentrations; (c) Phase diagram for the LCA boundary of surface-charge and salt-concentration with the 
same a on two planes; (d) Ionic densities when Dgff = 0.58 and \Alnm for systems of c^ = 0.1M from Panel 

(b). 
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tive pressure when the separation distance is comparable to Ed, for which the attractive 
strength decreases with the increase of cr. With the decrease of the bulk salt concentra¬ 
tion, the ionic density becomes smaller, reducing the entropic effect. This may reverse 
the attractive force at a certain concentration (shown in Panel (b)). We can see the at¬ 
traction occurs when the distance is at the scale of the Debye length for the two denser 
cases, which is consistent with the particle-based simulation results Il24ll . The pressure 
does not depend much on the asymmetry in surface charge densities, but the strength 
of the densities. By this observation, we will only study planes with the same a in the 
following. 

These above observations can be explained as that the pressure is resulted from the 
competition of the electrostatic repulsion and the entropically driven attraction IfTsH . and 
the LCA happens at the condition of smaller surface charge density for which the elec¬ 
trostatic repulsion is weak, or denser ionic density for which fhe entropically driven at¬ 
traction is strong. This is in agreement with Wang and Wang Il57ll although they consider 
the ions as point charges and can not reach systems with higher surface charges due fo 
the instability of the point-charge model IIMI . To study the behavior of higher surface 
charge density, we find the finite ion size is essential. The competition results in a phase 
boundary which separates the attractive region and the repulsive region in the (T—Ci, 
phase diagram, shown in Fig. |3](c). We see at the phase boundary, the relation is closes 
to Cfo ocu^, and accurate for small concentration up to C;, ~3M. The numerical results are 
shown for the bulk concentration up to ~6M, which illustrate the planes can be attractive 
for surface charge about O.le/ nm^, corresponding to a system in the mediated-coupling 
regime. The electrolyte with concentration as high as 6M is used in electrochemical ex¬ 
periments II33II . This square law can be understandable since the electrostatic interaction 
is proportional to at the weak-coupling regime, and the entropic contribution has the 
form of C;,log(c;,). 

To a closer look at the mechanism, in Fig. |3{d), we plot the ionic densities for systems 
of = O.IM and the surface charge densities = (72 = —O.Ole/nni^, corresponding to the 
solid circle line in Fig. |3]{b) with separation 0.58 and lAlnm. It shows clearly that when 
the planes are attractive, the ionic densities are much smaller than the bulk densities, 
leading to a significantly less pressure due to the ionic depletion. 

4.2 Asymmetric valences and nonuniform sizes 

The ionic size effect of ions is important for systems with nano-scale confinements. The 
inclusion of the ionic excluded volume and multivalent ions leads to a large change of 
the local correlation strength in self energy, and the consequence on the interface struc¬ 
ture and interaction is not well investigated. In this section, we perform numerical cal¬ 
culation for two groups of asymmefric sysfems. One group takes 1:1 electrolytes with 
variable ionic sizes; the other one takes uniform ionic sizes a± = O.lnm but varies the 
ionic valences. The surface charge densify and the bulk charge density remain constants, 
(7= —O.Ole/ nm^ and z+c+ = O.IM. The results are shown in Fig. HI 
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X {nm) X (nm) 


Figure 4: Pressure and charge distribution for asymmetric electrolytes with z+c^^O.lM between two planes of 
charge density <7= —0.Ole/nm^. (a) Pressure between two planes for 1:1 electrolytes of variable ionic sizes; (b) 
Pressure between two planes for asymmetric electrolytes and uniform ionic sizes; (c)(d) Total charge density 
distribution along the x-axis in electrolytes for (c) Dgff = 0.7nm and (d) Dgff = 1.9nm. 


The ionic size matters in two aspects: 1) the thickness of the Helmholtz layer where 
ions are inaccessible; and 2) the self energy. From Fig. HJa), we observe strong effects 
on pressure curves due to the non-uniformity of counterions and coions. For symmetric 
ionic sizes, the thickness of the Helmholtz layer is less important as we can see that the 
curves of a± = O.lnm and a± = OAnm are similar while the other two curves is totally 
different. For the effect of ionic size in the self energy, we should recall the formulation 
of the effective self energy, the difference between Uj and 


Ui — U^ 


M-|-Kfc/47r£oe 


(4.2) 


multiplied by the factor 2?e^/2. The leading contribution of u is about —K(r) plus the 
boundary contribution including the image charge. If the ion is not too close to the sur¬ 
face, u is in general a negative quantity, with—w < kj, due to the ion depletion. As a result. 
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the self-energy contribution of a larger ion is less than that of a smaller ion, and hence the 
a± = 0Anm curve behaves a weaker attraction than the a± = 0.2nm curve. If the ion size is 
very small which allows u > 0, then the image charge effect dominates in this area. The 
image charge effect becomes weaker as the ionic size gets larger, again ensuring a weaker 
attraction. 

For asymmetric ionic sizes, if counterions (cations) are bigger than coions (anions), the 
Helmholtz layer of counterions is thicker than coions'. Coions can be getting close to the 
charged surfaces, which effectively increases the density of the surface charge, thus the 
attractive force between interfaces is weakened or the interfaces become repulsive. This 
is shown by the curve of =0.2nm and =0.4nm. In contrast, if coions are bigger, the 
surface charge density is effectively attenuated, leading to a stronger attraction potential 
(see the curve of = 0Anm and = Q.2nm). In this case, the planes of short separation 
remain attractive. 

Let us look at the results in Fig. Stb), the effects of valence asymmetry. Since the self 
energy has quadratic relation with the valence, while ion-surface electrostatic force 

is roughly proportional to z„ the multivalent counterion receives very strong repulsion 
from the surface. If colons have smaller valence, which cause less repulsion, equiva¬ 
lently the surface charge density is increased, thus the attractive force is weakened. The 
results have shown this situation: the planes in 2:1 electrolytes have weaker attraction 
(becomes pure repulsion) than those in 1:1 and 1:2 electrolytes. The attractive distance 
also increases as the Debye length increases. 

The panels (c) and (d) of Fig. 0] present the total charge distribution along the nor¬ 
mal direction of the interfaces for two separations Dgff = 0.7nm and 1.9nm from panel 
(b). It can be observed that the 2:-l electrolyte does show negative value near interfaces, 
validating the preceding analysis that the high-valence counterions increase the effective 
surface charge density. It should be mentioned that the curve for the 1:2 electrolyte in 
panel (d) is concave, different from the other three curves. Near the middle point, the 
total charge reverses the sign, illustrating a charge inversion, another important elec¬ 
trostatic many-body phenomenon IIT0ll26ll . Again, it is owing to the strong repulsion of 
multivalent anions (counterions) by their image charge. 

4.3 Dielectric variation effect 

We calculate the pressure between planes by varying the dielectric constant of the me¬ 
dia outside the surfaces, eg, from 2.5 to 160. The results given in Fig. |5la) illustrate a 
monotonic increase of the minimum pressure with eg. Compared to the homogeneous 
permittivity eg = 80, the low dielectric constant of the boundary media greatly enhances 
the LCA. Interestingly, the planes become always repulsive when eg is much larger than 
the water dielectric constant, e.g., eg = 160. This shows a counter-intuitive phenomenon 
that the conducting limit of the planes tends to be repulsive due to the entropic forces. 

In Fig. ISlJb), we show the comparison with a variable permittivity in water. Since wa¬ 
ter molecule is polarizable, the dielectric permittivity of the water medium depends on 
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Figure 5; Effects on dielectric permittivity on the pressure between two planes. The 1:1 electrolytes with 
a±=0.2nm (a) Water dielectric constant ew = 80 and varying boundary dielectric constant eg. Surface charge 
tr = —O.Ole/wm^, bulk concentration 0.05M; (b) Boundary dielectric constant eg =2.5 and field-dependent 
dielectric permittivity for water. Surface charge cr=—0.03e/nm^, concentration 0.5M. 


the strength of the electrostatic field due to the orientation of dipole molecules against the 
field. A more ordered orientation leads to smaller permittivity, which can be described 
by the Langevin model II23II . 

e(|VO|) = l+3(£w-l)>C(/3po|VO|)//lpo|VO|, (4.3) 

where C{y) = coth(y) — l/y is the Langevin function, and po is the permanent dipole 
moment. We take po=4.8D and ew=80. With the weak surface charge density, the results 
show the electric field near the surfaces slightly increases the interface-interface pressure 
which is due to a stronger electrostatic interaction in lower dielectric permittivity. 


5 Conclusions 

In summary, we study the phenomenon of the LCA between charged planes from a 
recently-developed self-consistent field model. We investigate the influence of different 
parameters of the electrical double layer on the many-body phenomenon. It is shown the 
depletion-induced LCA depends on the surface charge density and the bulk salt concen¬ 
tration, and the dielectric mismatch significantly changes the pressure strength between 
two planes. We find the mechanism of inducing LCA is mostly the depletion effect be¬ 
tween the charged surfaces. It is more likely to observe LCA for low surface charge, high 
salt concentration and the attraction can be greatly enhanced by small permittivity of the 
dielectric media outside the electrolyte. However, as the distance between charged sur¬ 
faces which have LCA is comparable with Debye length, high salt concentration might 
lead to very narrow separation. 
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